Supplement A: Estimating Travel Cost
1 Overview
Goal: estimate cost-distance (in time) from each ranch to Ciudad Constitucion, La Paz, and local springs.
Data: a Digital Elevation Model (DEM) along with point locations for ranche clusters, springs, cities.
Method: Campbell’s hiking function (Campbell et al. 2019).
There are two watersheds in the project area. Each one has a single road that connects it to Highway 1 and thus to the markets in Ciudad Constitucion to the north and the capitol La Paz to the south. The “road” that connects the watersheds to each other is not maintained and thus rarely used. As a consequence, nearly all of the variance in travel time to cities is found within each watershed. For that reason, we estimate the time it takes to leave each watershed (“leaving” here means hitting the edge of the project area), rather than the true time it takes to get to each city. We interpret this as the relative difference in time it takes to get to Ciudad Constitucion from the northern watershed and the relative difference in time it takes to get to La Paz from the southern watershed. To get to La Paz from the northern watershed, we simply add an arbitrary half hour to the total time to get out of that watershed. The same would presumably be the case to get to Ciudad Constitucion from the southern watershed, though that value does not factor into our analysis.
For springs, Dr. MacFarlan’s ethnographic data suggest that individuals residing within the same ranch cluster tend to rely on one or two springs at most. Thus, we estimate road distance to the two nearest springs and take the average.
All spatial data are projected to the Mexico ITRF92 / UTM zone 12N CRS (EPSG:4485).
2 R Preamble
library(dplyr)
library(ggnewscale)
library(ggplot2)
library(here)
library(igraph)
library(sf)
library(terra)
library(viridis)Specify consistent plot theme here:
Code
# this is designed to position the legend
# at the top left corner of the plot area
theme_set(theme_bw(12))
theme_update(
axis.text = element_text(size = rel(0.5), color = "gray45"),
axis.text.y = element_text(angle = 90, hjust = 0.5),
axis.ticks = element_line(color = "gray45", linewidth = 0.1),
axis.ticks.length = grid::unit(0.1, "cm"),
axis.title = element_blank(),
legend.background = element_rect(fill = "transparent", color = "transparent"),
legend.key.height = grid::unit(0.4, "cm"),
legend.key.width = grid::unit(0.3, "cm"),
legend.margin = margin(l = 0),
legend.justification = c("left", "top"),
legend.spacing.x = grid::unit(0.05, "cm"),
legend.spacing.y = grid::unit(0.05, "cm"),
legend.text = element_text(vjust = 0.7, margin = margin()),
legend.title = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(size = rel(1.1), face = "bold")
)
point_colors <- c("#FDE74C", "#EE7733", "#0077BB")
pretty_breaks <- function(x, .axis = "x", dx = 2500, n = 4){
mn <- paste0(.axis, "min")
mx <- paste0(.axis, "max")
brks <- seq(x[[mn]] + dx, x[[mx]] - dx, length.out = n)
round(brks, digits = -3)
}3 Data
3.1 Geopackage
Specify path to geopackage database holding all spatial vector data.
gpkg <- here("data", "choyero.gpkg")3.2 Project Window
bcs <- read_sf(gpkg, layer = "baja")
window <- read_sf(gpkg, layer = "window")
roads <- read_sf(gpkg, layer = "roads")
watersheds <- read_sf(gpkg, layer = "watersheds")3.3 Elevation (DEM)
dem <- rast(here("data", "dem_30m.tif"))3.4 O-D Points
# ORIGIN ---
clusters <- gpkg |>
read_sf(layer = "clusters") |>
mutate(variable = "cluster")
# DESTINATION --- to city, terminus at edge of window
terminus <- gpkg |>
read_sf(layer = "terminus") |>
mutate(variable = "terminus")
springs <- gpkg |>
read_sf(layer = "springs") |>
select(id, name) |>
mutate(variable = "spring")Code
od_points <- clusters |>
rename("name" = cluster) |>
mutate(name = as.character(name)) |>
select(name, variable) |>
bind_rows(terminus, springs) |>
mutate(variable = ifelse(variable == "terminus", "city", variable))
bb8 <- st_bbox(window)
dx <- bb8[["xmax"]] - bb8[["xmin"]]
dy <- bb8[["ymax"]] - bb8[["ymin"]]
# add hillshade relief to visualize elevation
# https://dominicroye.github.io/en/2022/hillshade-effects/
hillshade <- shade(
slope = terrain(dem, "slope", unit = "radians"),
aspect = terrain(dem, "aspect", unit = "radians"),
angle = 45,
direction = seq(45, 360, by = 45),
normalize = TRUE
)
hillshade <- mask(sum(hillshade), vect(st_union(watersheds)))
names(hillshade) <- "shade"
gg <- ggplot() +
geom_raster(
data = as.data.frame(hillshade, xy = TRUE),
aes(x, y, fill = shade)
) +
scale_fill_distiller(
palette = "Greys",
na.value = "transparent",
guide = "none"
) +
new_scale_fill() +
geom_raster(
data = as.data.frame(dem, xy = TRUE),
aes(x, y, fill = elevation),
alpha = 0.7
) +
scale_fill_distiller(
palette = "Greys",
na.value = "transparent",
guide="none"
) +
new_scale_fill() +
geom_sf(
data = watersheds,
fill = "transparent",
color = alpha("darkblue", 0.35),
linewidth = 0.2
) +
geom_sf(
data = roads,
color = "black",
linewidth = 0.3
) +
geom_sf(
data = od_points,
aes(shape = variable, fill = variable, color = variable),
stroke = 0.3,
size = 1.8
) +
scale_color_manual(values = c("grey15", "white", "white")) +
scale_fill_manual(values = alpha(point_colors, 0.75)) +
scale_shape_manual(values = c(23, 22, 21)) +
coord_sf(datum = 4485) +
scale_x_continuous(
breaks = pretty_breaks(bb8, "x"),
expand = expansion(add = 1000)
) +
scale_y_continuous(
breaks = pretty_breaks(bb8, "y", n = 3),
expand = expansion(add = 1000)
) +
theme(legend.position = c(0.03, 0.98))
ggsave(
plot = gg,
filename = here("figures", "od-points.png"),
dpi = 800,
width = 3.5,
height = 3.5 * (dy/dx)
)
4 Path analysis
To estimate least-cost paths, we first derive slope estimates from the DEM then calculate the hiking speed at which individuals can traverse any given cell using Campbell’s hiking function. After that, we decrease the cost in units of time to traverse grid cells crossed by roads and other paths using weights that approximate travel times by car or horse. This is equivalent to increasing the speed one travels in those cells. Below are the specific weights we use.
| Type | Weight | Approx. Speed |
|---|---|---|
| primary | 0.2 | 25 mph (41.5 kmh) |
| secondary | 0.5 | 10 mph (16.2 kmh) |
| tertiary | 0.8 | 07 mph (11.0 kmh) |
The necessary functions are defined in the folded code chunk below.
Code
#' Cost surface
#'
#' Create a `terrain` representing the cost of travel.
#'
#' The `SpatRaster` must have a projected CRS, with both distance units and elevation
#' values in meters.
#'
#' @param x a `SpatRaster` with elevation values.
#' @param hf a character string specifying the preferred hiking function. Options
#' include `tobler` and `campbell` (default).
#' @param max_slope a numeric value specifying the maximum allowed hiking slope
#' _in degrees_ (default is 45). Suppose you ask: Would a person really hike
#' _that_? This is the shallowest slope where the answer to that question is No.
#' @param neighbors numeric, one of 4, 8, or 16 (default) specifying the number of
#' adjacent cells between which to calculate movement cost. These are passed on
#' to `terra::adjacent()` and represent standard neighborhoods (rook, queen, and
#' knight + queen).
#' @param ... Arguments passed on to the hiking function specified by `hf`.
#'
#' @return a `terrain` representing cost of travel.
hf_terrain <- function(x,
hf = "campbell",
max_slope = 45,
neighbors = 8,
...) {
# find adjacent cells in the raster
# make sure to ignore all na values
n_cells <- terra::ncell(x)
na_cells <- which(is.na(terra::values(x)))
cells <- setdiff(1:n_cells, na_cells)
adj <- terra::adjacent(
x,
cells = cells,
directions = neighbors,
pairs = TRUE
)
adj <- adj[!adj[, 2] %in% na_cells, ]
### slope ###
# change in y between adjacent pixels
heights <- terra::values(x)[, 1]
rise <- (heights[adj[, 2]] - heights[adj[, 1]])
# change in x between adjacent pixels
run <- terra::distance(
terra::xyFromCell(x, adj[, 1]),
terra::xyFromCell(x, adj[, 2]),
lonlat = FALSE,
pairwise = TRUE
)
slope <- rise/run
### speed ###
speed_fun <- switch(
hf,
"campbell" = campbell,
"tobler" = tobler,
stop(paste0("hiker does not implement the ",
hf,
" hiking function."))
)
speed <- speed_fun(slope, ...)
### conductance ###
conductance <- speed/run
# what is a realistic slope people would try to hike?
max_slope <- tan(max_slope * pi / 180) # degrees -> radians -> rise-over-run
i <- which(abs(slope) >= max_slope)
conductance[i] <- 0
# build matrix
cm <- Matrix::spMatrix(nrow = n_cells, ncol = n_cells)
cm[adj] <- conductance
# include some spatial and other information to make
# manipulating the network faster and easier
# when passed to other functions
# storing bounding box as vector to keep the size down
bb8 <- terra::ext(x)
bb8 <- c(bb8$xmin, bb8$xmax, bb8$ymin, bb8$ymax)
xcrs <- terra::crs(x, describe = TRUE)
xcrs <- paste0(xcrs$authority, ":", xcrs$code)
tmin <- 1/max(conductance)
tmax <- 1/min(conductance[conductance > 0])
terrain <-
list(
"conductance" = cm,
"neighbors" = neighbors,
"nrow" = terra::nrow(x),
"ncol" = terra::ncol(x),
"bb8" = bb8,
"crs" = xcrs,
"range" = c("min" = tmin, "max" = tmax)
)
class(terrain) <- "terrain"
return(terrain)
}
campbell <- function(x, decile = 50L) {
if (length(decile) != 1) {
stop("campbell accepts only one decile at a time.")
}
# campbell function wants degrees, not rise-over-run as in tobler
# but, we're working with rr, so remember to convert radians to degrees with 180/pi!
slope <- atan(x) * 180 / pi
# values provided in Campbell 2019, Supplement, Table S1
# for simplicity, I excluded all but the deciles
deciles <-
data.frame(
decile = c(10L, 20L, 30L, 40L, 50L, 60L, 70L, 80L, 90L),
a = c(-1.568, -1.71, -1.858, -1.958, -2.171, -2.459, -2.823, -3.371, -3.06),
b = c(13.328, 10.154, 8.412, 8.96, 10.064, 11.311, 12.784, 15.395, 16.653),
c = c(38.892, 36.905, 39.994, 50.34, 63.66, 79.287, 98.697, 134.409, 138.875),
d = c( 0.404, 0.557, 0.645, 0.649, 0.628, 0.599, 0.566, 0.443, 0.823),
e = c(-0.003, -0.004, -0.004, -0.005, -0.005, -0.005, -0.005, -0.005, -0.0139)
)
ind <- which(deciles$decile == as.integer(decile))
deciles <- deciles[ind, ]
a <- deciles$a
b <- deciles$b
c <- deciles$c
d <- deciles$d
e <- deciles$e
# lorentz distribution
lorentz <- (1 / ((pi * b) * (1 + ((slope - a) / b)^2)))
# modified lorentz
(c * lorentz) + d + (e * slope)
}
#' Channel of reduced cost
#'
#' Reduces cost of travel along paths that cross through channels or
#' other aspects of a terrain that expedite travel.
#'
#' @param x a cost `terrain` generated by `hf_terrain()`.
#' @param channel an `sf` object.
#' @param .m a numeric between 0 and 1 applied as a weight to the cost of travel.
#'
#' @return a `terrain` representing cost of travel.
hf_channel <- function(x, channel, .m = 1) {
rr <- terra::rast(
nrow = x$nrow,
ncol = x$ncol,
extent = terra::ext(x$bb8),
crs = x$crs
)
channel <- terra::vect(channel)
cells <- terra::cells(rr, channel)[, 2]
adj <- cbind(
"from" = as.integer(x$conductance@i + 1),
"to" = as.integer(x$conductance@j + 1)
)
# find incident edges, from or to cells intersecting barrier
i <- which(adj[, 1] %in% cells | adj[, 2] %in% cells)
adj <- adj[i, ]
# want to reduce time, which means dividing conductance by .m
x$conductance[adj] <- (x$conductance[adj] / .m)
return(x)
}primary_roads <- roads |> filter(level == "P") |> st_buffer(dist = 30)
secondary_roads <- roads |> filter(level == "S") |> st_buffer(dist = 30)
tertiary_roads <- roads |> filter(level == "T") |> st_buffer(dist = 30)
terrain <- dem |>
hf_terrain(max_slope = 35, decile = 30) |>
hf_channel(primary_roads, .m = 0.2) |>
hf_channel(secondary_roads, .m = 0.5) |>
hf_channel(tertiary_roads, .m = 0.8)
remove(primary_roads, secondary_roads, tertiary_roads)Code
#' Terrain to raster
#'
#' Converts a cost `terrain` into a raster representing
#' travel cost, specifically the average cost of moving from/to each
#' cell. Mostly useful for visualizing with `terra::plot()`
#' as a sanity check.
#'
#' @param x a cost `terrain` generated by `hf_terrain()`.
#'
#' @return a `SpatRaster` with travel cost values.
hf_rasterize <- function(x) {
rr <- terra::rast(
nrow = x$nrow,
ncol = x$ncol,
extent = terra::ext(x$bb8),
crs = x$crs
)
# this part is basically what Jacob van Etten does with
# gdistance::raster(TransitionLayer), though that offers more options
# here I am just averaging the vertical and horizontal sums
total_row <- Matrix::rowSums(x$conductance)
total_col <- Matrix::colSums(x$conductance)
#logical sparse matrix, identifies number of adjacent cells
lm <- methods::as(x$conductance, "lMatrix")
n_row <- Matrix::rowSums(lm)
n_col <- Matrix::colSums(lm)
# average cost of moving from/to each cell
# values <- total_cost / n_adjacent
v_row <- total_row / n_row
v_col <- total_col / n_col
values <- (v_row + v_col)/2
terra::setValues(rr, (1/values))
}
gg <- ggplot() +
geom_raster(
data = as.data.frame(hf_rasterize(terrain), xy = TRUE),
aes(x, y, fill = lyr.1)
) +
scale_fill_viridis(
name = "Seconds",
breaks = c(10, 30, 50)
) +
coord_sf(datum = 4485) +
scale_x_continuous(
breaks = pretty_breaks(bb8, "x"),
expand = expansion(add = 1000)
) +
scale_y_continuous(
breaks = pretty_breaks(bb8, "y"),
expand = expansion(add = 1000)
) +
guides(
fill = guide_colorbar(
barwidth = 3.8,
barheight = 0.6,
raster = FALSE,
frame.linewidth = 0.2,
ticks.linewidth = 0.7,
title.position = "top"
)
) +
theme(
legend.position = c(0.03, 0.98),
legend.direction = "horizontal",
legend.spacing.x = grid::unit(0.05, "cm"),
legend.text = element_text(size = rel(0.5), margin = margin(l = 2)),
legend.title = element_text(
size = rel(0.8),
margin = margin(b = 2, t = 2),
hjust = 0
)
)
ggsave(
plot = gg,
filename = here("figures", "cost-surface.png"),
dpi = 600,
width = 3.5,
height = 3.5 * (dy/dx)
)
4.1 To springs
Function to calculate paths and travel times is in the folded code chunk below.
Code
#' Shortest path between points
#'
#' Calculate the "shortest" or least-cost paths between multiple start and end points.
#' For each from-point, the shortest path to every to-point is calculated.
#'
#' @param x a cost `terrain` generated by `hf_terrain()`.
#' @param from an `sf` specifying start locations, must be POINT or MULTIPOINT.
#' @param to an `sf` specifying end locations, must be POINT or MULTIPOINT.
#' @param add_cost logical, whether to include a column of travel costs.
#'
#' @return An `sf` object with least-cost `LINE` features and two columns: from and to.
#' Values in the from- and to-columns are `1:nrow(from)` and `1:nrow(to)`, respectively.
#' If `add_cost = TRUE`, the `sf` object also includes a travel cost column.
hf_hike <- function(x, from, to, add_cost = FALSE) {
from_xy <- sf::st_coordinates(from)[, 1:2, drop = FALSE]
to_xy <- sf::st_coordinates(to)[, 1:2, drop = FALSE]
rr <- terra::rast(
nrow = x$nrow,
ncol = x$ncol,
extent = terra::ext(x$bb8),
crs = x$crs
)
from_cells <- terra::cellFromXY(rr, from_xy)
to_cells <- terra::cellFromXY(rr, to_xy)
graph <- igraph::graph_from_adjacency_matrix(
x$conductance,
mode = "directed",
weighted = TRUE
)
# invert conductance to get travel cost
igraph::E(graph)$weight <- (1 / igraph::E(graph)$weight)
# return a list of lists of vectors of vertices on paths
shorties <-
lapply(
from_cells,
function(z) {
res <-
igraph::shortest_paths(
graph,
from = z,
to = to_cells,
mode = "out"
)
res$vpath
})
# collapse (list of lists of vectors) to (list of vectors)
shorties <- unlist(shorties, recursive = FALSE)
# now we can just make a bunch of linestrings out of them
line_list <-
lapply(
shorties,
function(z) {
cells <- as.vector(z)
# to handle case where from-vertex == to-vertex
if (length(cells) == 1) cells <- rep(cells, 2)
xy <- terra::xyFromCell(rr, cells)
sf::st_linestring(xy)
}
)
# and then an sf
sf_col <- sf::st_sfc(line_list, crs = x$crs)
short_paths <- sf::st_sf(geometry = sf_col)
# add indices
short_paths <- transform(
short_paths,
from = rep(1:nrow(from_xy), each = nrow(to_xy)),
to = rep(1:nrow(to_xy), times = nrow(from_xy))
)
if (add_cost) {
cost <- igraph::distances(graph, from_cells, to_cells, mode = "out")
short_paths <- transform(short_paths, cost = c(t(cost)))
}
return(short_paths)
}to_springs <- hf_hike(
terrain,
from = clusters,
to = springs,
add_cost = TRUE
)
to_springs <-
to_springs |>
group_by(from) |>
slice_min(cost, n = 2) |> # for each cluster, get two lowest cost paths
mutate(cost = (cost/3600)) |> # convert to hours
summarize(
cluster = unique(from),
spring = paste(to, collapse = ","),
cost = mean(cost)
) |>
select(cluster, spring, cost)4.2 To cities
Please note that travel time is not strictly to cities, but to the edge of the project area. As noted above, for the northern watershed, the time to La Paz equals the time it takes to leave the watershed plus half an hour.
north <- clusters |>
st_filter(
st_union(filter(watersheds, id %in% c(3294, 3296, 3299)))
)
south <- clusters |>
st_filter(
filter(watersheds, id == 3326)
)
### NORTH ###
out_north <- hf_hike(
terrain,
from = north,
to = terminus[1, ],
add_cost = TRUE
) |>
mutate(
cluster = north$cluster,
city = "Ciudad Constitucion",
cost = (cost/3600),
watershed = "north"
)
### SOUTH ###
out_south <- hf_hike(
terrain,
from = south,
to = terminus[2, ],
add_cost = TRUE
) |>
mutate(
cluster = south$cluster,
city = "La Paz",
cost = (cost/3600),
watershed = "south"
)
to_cities <-
bind_rows(out_north, out_south) |>
mutate(time_to_paz = ifelse(watershed == "north", cost + 0.5, cost)) |>
select(cluster, city, cost, time_to_paz)
remove(north, south, out_north, out_south)4.3 LCP Map
Code
od_points <- bind_rows(
od_points |> mutate(destination = ifelse(variable == "spring", "To springs", "To cities")),
od_points |> filter(variable == "cluster") |> mutate(destination = "To springs")
)
paths <- bind_rows(
to_springs |> mutate(destination = "To springs"),
to_cities |> mutate(destination = "To cities")
)
gg <- ggplot() +
geom_raster(
data = as.data.frame(hillshade, xy = TRUE),
aes(x, y, fill = shade)
) +
scale_fill_distiller(
palette = "Greys",
na.value = "transparent",
guide = "none"
) +
new_scale_fill() +
geom_raster(
data = as.data.frame(dem, xy = TRUE),
aes(x, y, fill = elevation),
alpha = 0.7
) +
scale_fill_distiller(
palette = "Greys",
na.value = "transparent",
guide="none"
) +
new_scale_fill() +
geom_sf(
data = watersheds,
fill = "transparent",
color = alpha("darkblue", 0.35),
linewidth = 0.2
) +
geom_sf(
data = roads,
color = "black",
linewidth = 0.3
) +
geom_sf(
data = paths |> st_buffer(125) |> group_by(destination) |> summarize(),
fill = "#C7FFED",
color = "#00291D",
linewidth = 0.2
) +
geom_sf(
data = od_points,
aes(shape = variable, fill = variable, color = variable),
stroke = 0.3,
size = ifelse(od_points$variable == "city", 2.2, 1.5)
) +
facet_wrap(vars(destination)) +
scale_color_manual(values = c("grey15", "white", "white")) +
scale_fill_manual(values = alpha(c(point_colors, 0.75))) +
scale_shape_manual(values = c(23, 22, 21)) +
coord_sf(datum = 4485) +
scale_x_continuous(
breaks = pretty_breaks(bb8, "x"),
expand = expansion(add = 1000)
) +
scale_y_continuous(
breaks = pretty_breaks(bb8, "y", n = 3),
expand = expansion(add = 1000)
) +
theme(
legend.position = c(0.01, 0.98),
strip.background = element_blank(),
strip.text = element_text(size = rel(1.1), hjust = 0)
)
ggsave(
plot = gg,
filename = here("figures", "least-cost-paths.png"),
dpi = 600,
width = 5.75,
height = 5.75 * (dy/(2*dx)) + 0.47
)
4.4 Save
write_sf(
to_springs,
dsn = gpkg,
layer = "path_to_springs"
)
write_sf(
to_cities,
dsn = gpkg,
layer = "path_to_cities"
)5 References
Campbell, Michael J., Philip E. Dennison, Bret W. Butler, and Wesley G. Page. 2019. “Using Crowdsourced Fitness Tracker Data to Model the Relationship Between Slope and Travel Rates.” Applied Geography 106: 93–107. https://doi.org/https://doi.org/10.1016/j.apgeog.2019.03.008.
6 Session Info
Code
# save the session info as an object
pkg_sesh <- sessioninfo::session_info(pkgs = "attached")
# inject the quarto info
pkg_sesh$platform$quarto <- paste(
quarto::quarto_version(),
"@",
quarto::quarto_path()
)
# print it out
pkg_sesh─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.3.1 (2023-06-16 ucrt)
os Windows 11 x64 (build 22621)
system x86_64, mingw32
ui RTerm
language (EN)
collate English_United States.utf8
ctype English_United States.utf8
tz America/Denver
date 2024-01-23
pandoc 3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
quarto 1.4.545 @ C:\\PROGRA~1\\Quarto\\bin\\quarto.exe
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.3.1)
ggnewscale * 0.4.9 2023-05-25 [1] CRAN (R 4.3.1)
ggplot2 * 3.4.4 2023-10-12 [1] CRAN (R 4.3.2)
here * 1.0.1 2020-12-13 [1] CRAN (R 4.3.1)
igraph * 1.6.0 2023-12-11 [1] CRAN (R 4.3.2)
sf * 1.0-15 2023-12-18 [1] CRAN (R 4.3.2)
terra * 1.7-65 2023-12-15 [1] CRAN (R 4.3.1)
viridis * 0.6.4 2023-07-22 [1] CRAN (R 4.3.1)
viridisLite * 0.4.2 2023-05-02 [1] CRAN (R 4.3.1)
[1] C:/Users/kenne/AppData/Local/R/win-library/4.3
[2] C:/Program Files/R/R-4.3.1/library
──────────────────────────────────────────────────────────────────────────────